Multi-product expansion for Nonlinear 
Differential Equations 



Jiirgen Geiser 
juergen.geiser@uni-greifswald.de 

Abstract. In this paper we discuss the extention of MPE methods to 
nonhnear differential equations. We concentrate on nonhnear systems of 
differential equations and generalize the recent MPE method, see |16) . 

Keyword Suzuki's method, Multi-product expansion, exponential splitting, 
non-autonomous systems, nonlinear differential equations. 

AMS subject classifications. 65M15, 65L05, 65M71. 
1 Introduction 

In this paper we concentrate on approximation to the solution of the nonlinear 
evolution equation, e.g. nonlinear evolution equation, 

dtu^ F{u{t)) = A{u{t)) + B{u{t)), u{0) = uo, (1) 

with the unbounded operators A : C X X and B : D{B) C X ^ X. 

We have further F{v) = A{v) + B{v), v G D{A) D D{B). 

We assume to have suitable chosen subspaces of the underlying Banach space 
(X, II • llx) such that D{F) = D{A) n D{B) ^ 

For such equations, we concentrate on comparing the higher order meth- 
ods to Suzuki's schemes. Here the Suzuki's methods apply factorized symplectic 
algorithms with forward derivatives, see [TT], [T^. 

The exact solution of the evolution problem ^ is given as: 

u{t) = £F{t, u(0)), < t < T, (2) 

with the evolution operator Bp depending on the actual time t and the initial 
value m(0). 

We apply a formal notation given as: 

u(t) = exp(iDF)M(0), < t < T, (3) 

Here the evolution operator exp(il?i?) and the Lie-derivative Dp associated 
with F are given as: 

%y.^{tDF)Gv = GiSpit, v)), 0<t<T, DpGv = G'{V)F{v) (4) 

for any unbounded nonlinear operator G : D{G) C X X with Frechet deriva- 
tive G'. 
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2 Nonlinear Exponential operator splitting methods 

In the course of devising numerical algorithms for solving the prototype nonlinear 
differential equations 

dtu = A{u) + B{u), u(0) = uo, (5) 

where A and B are non-commuting operators, 

Strang [38] proposed two second-order algorithms corresponding to approxi- 
mating 

T{At) = e'^*(^-^+^«) (6) 

either as 

S{At) = i (e^*^^e^*^- + e^*^«e^*^^) (7) 

or as 

SAB{At) - c^^tmDn^AtD,^iAt/2)Ds^ (g) 

Following up on Strang's work, Burstein and Mirin[S] suggested that Strang's 
approximations can be generalized to higher orders in the form of a multi-product 
expansion (MPE), 

k i 

and gave two third-order approximations 



3 Y 2 / 3 

where Sba is just Sab with A -k^ B interchanged, and 



(10) 



BAB{At) = ^QiAt/S)DA^i2At/3)DB^{2At/3)DA^{At/3)DB _ l^AtDA^AtDs _ (jj^ 

8 8 

They credited J. Dunn for finding the decomposition D{At) and noted that the 
weights Cfe are no longer positive beyond second order. Thus the stability of the 
entire algorithm can no longer be inferred from the stability of each component 
product. 



3 Nonlinear Multi-product decomposition 

The multi-product decomposition ([9]) is obviously more complicated than the 
single product splitting. 

By the way after Burstein and Mirin, Sheng[57[ proved their observation that 
beyond second-order, Uki, bki and Ck cannot all be positive. Some possibilities 
are the idea of complex coefficients to obtain a higher order scheme. 

For general applications, including solving time-irreversible problems, one 
must have ati and bki positive. 
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Here we show, that we obtain a extrapolation scheme, that can overcome 
such problems. 

Therefore every single product in (|9]) can at most be second-order f3 7|39] . 
But such a product is easy to construct, because every left-right symmetric 
single product is second-order. Let Ts{h) be such a product with "^fej ~ 1 
and Xj bki — 1, then Tsih) is time-symmetric by construction, 

Ts{~h)Ts{h) = 1, (12) 

implying that it has only odd powers of h 

Ts{At) = eMMDA + Db) + At^Ej, + At'^E^ + • • • ) (13) 

and therefore correct to second-order. (The error terms Ei are nested commu- 
tators of Da and Db depending on the specific form of Ts.) This immediately 
suggests that the fcth power of 7s at step size h/k must have the form 

TsiAt/k) = exp{At{DA + Db) + k'^At^E^i + k'^At^E^ + ■■■), (14) 

and can serve as a basis for the multi-production expansion ([9]). The simplest 
such symmetric product is 

T2{h) = SABih) or T2{h) =^ SBA{h). (15) 

If one naively assumes that 

r2{h) = Q^^iD^+Ds) + Ch^ + Dh'^ + ■■■ , (16) 

then a Richardson extrapolation would only give 

— ^ [k'Tiih/k) ~ T2{h)\ = e^*(^^+^^) + 0(/i4), (17) 

a third-order [35] algorithm. However, because the error structure of T2{h/k) is 
actually given by (fT4|) . one has 

T^{h/k)^e'^'^''-+''-Kk-^h''E^+]^k-^h%DA+DB)E^+E^{DA+DB)]+0{h'>), 

(18) 

and both the third and fourth order errors can be eliminated simultaneously, 
yielding a fourth-order algorithm. Similarly, the leading 2n + l and 2n-\-2 order 
errors are multiplied by fc^^" and can be eliminated at the same time. Thus for 
a given set of n whole numbers {ki\ one can have a 2rith-order approximation 

^At(D.+Ds) ^ J2 C.r/-' + 0(/l^"+l). (19) 

provided that Ci satisfy the simple Vandermonde equation: 
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Surprising, this equation has closed form solutions [13] for all n 

The natural sequence {fc,} = {1, 2, 3 ... n} produces a 2nth-order algorithm with 
the minimum n{n + l)/2 evaluations of T2{h). For orders four to ten, one has 
explicitly: 

%iAt) = -lT2{At) + ^T^(^) (22) 



3 ^ V 2 

1 16^2 , 81 3 fAt\ 



%iAt) = -TMt) - -r,^ ^- J + -r/ ^- j (23) 

16384^4 /zit\ 390625 5 ( At\ 



2835 J ^ 72576"'^2 J " ^^^^ 

As shown in Ref.[13). Ti{h) reproduces Nystrom's fourth-order algorithm with 
three force-evaluations and Teih) yielded a new sixth-order Nystrom type algo- 
rithm with five force-evaluations. 

3.1 Error Analysis of the Nonlinear Mult ipro duct Expansion 

The error analysis is based on the following ideas: 

— Definition of a new operator (linear operator) 

— Definition of the derivatives 

— Derivation of the splitting errors 

The nonlinear equation is given as: 

dtu^ F{u{t)) = A{u{t)) + B{u{t)), u{0) = uo, (26) 

We associate a new operator F which is a linear Lie operator: 

Definition 1. For the given operator F we associate a new operator, denoted 
F. 

This Lie-operator acts on the space of the differentiable operators of the type: 
S^S , 

and maps each operator G into the new operator F{G), such that for any 
element c E S : 

(F(G))(c) = (G'(c)oF)(c), (27) 
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Here the derivatives are given as: 

Definition 2. The k-th power of the Lie-operator F applied to some operator 
G can he expressed as the k-th derivative of G, that is the relation 

{fHG)){c) = (28) 

is valid for all fc = 1, 2, . . .. 

Let us splitt the operator F into the sum Fi + F2 
The sphtting error for the A-B sphtting is given as: 

ErrA-B = exp(TA + F2){I) - (exp(rFi) exp(rF2))(/)(c), (29) 

if we set in the commutator we obtain for Fi and F2 

ErrA-BX2 = riF^ic) o Fi)(c) - Fi{c) o F2){c)) + 0{t^), (30) 

The Strang Sphtting is given as: 

Errstrang = exp(rFi + A ) (/) (31) 
-(exp(T/2Fi) exp(TF2) exp(T/2A))(/)(c), 

if we set in the commutator we obtain for Fi and F2 

ErrstrangA.2 = ^^'([^2, [-^2 , ^1 ] ] ( c) - 2[Fi, [F,,F2]]{c)) + 0{t^) . (32) 

with the commutator: [Fi,F2]{c) = (F^(c) o Fi)(c) - {Fi{c) o i^2)(c). 

We derive the nonlinear MPE based on the definition of the hnearized Lie- 
operator. 

The derivation of the MPE method is given as: 

r2^{h/k) = e'''^''^+''^^+k-^h^E3+^k-^h''[iDA+DB)E3+E3{DA+DB)]+Oih^), 

(33) 

and both the third and fourth order errors can be ehminated simuhaneously, 
while E3 and E4 are given as commutators of the derivatives of A and B, see 

m- 

4 Nonlinear Operator splitting methods 

In the literature, there are various types of splitting methods. We mainly con- 
sider the following operators splitting schemes in this study: 



6 



1. Sequential operator splitting: A-B splitting 

^'^*^^'' =A{c*{t)) witht e and c*(r)=c" (34) 



at ' ' - L > J y J sp 

B{c**{t)) withi e and c**{t") = c*{e+^), (35) 



dt 

for n = 0, 1, iV — 1 whereby c" = cq is given from (??). The approximated 
split solution at the point t = is defined as Cg+^ = c**(i"+^). 
2. Strang-Marchuk operator splitting : A-B- A splitting 



= A(c*(t)) withi e and c*(r) = c^p (36) 

^-^^^ =B(c**(t)) withi e [i",t"+V2] and c"(r) =c*(r+V2) , (37) 
= A{c*{t)) with t G and c***(f"+i/2) = c**{t''+^^S) 

where t"+'^/^ = t"' + t/2, t is the local time step. The approximated split solu- 
tion at the point t = is defined as c^+^ = c***(r+i). 

3. Iterative splitting with respect to one operator 



dr(t) 

-^=A{ci{t)) + B(Q_i(i)), with din = c",i = 1,2,..., m (39) 
4. Iterative splitting with respect to alternating operators 



^""'^^^ = A{c,{t)) + B{ci-i{t)), with c,(t") = c", (40) 



dt 

i = l,2,...,j , 

= A{c,-i{t)) + B{c,{t)), with c,+i(r) = c", (41) 

= J + 1, j + 2,. . . ,TO , 

Here, co(t") = c" , c_i = and c" is the known split approximation at the 
time level t = t". The split approximation at the time-level t = is defined 
asc"+i =C2™+i(i"+i). 



4.1 Error Analysis of the Classical Splitting Schemes 

We associate a new operator F which is a linear Lie operator: 



7 



Definition 3. For the given operator F we associate a new operator, denoted 
F. 

This Lie-operator acts on the space of the differentiable operators of the type: 
S^S , 

and maps each operator G into the new operator F[G), such that for any 
element c ^ S: 

(F(G))(c) = (G'(c)oF)(c), (42) 
Here the derivatives are given as: 

Definition 4. The k-th power of the Lie-operator F applied to some operator 
G can he expressed as the k-th derivative of G, that is the relation 

iF\G)){c) = (43) 

is valid for all fc = 1, 2, . . .. 

Let us splitt the operator F into the sum Fi + F2 
The splitting error for the A-B splitting is given as: 

ErrA-B = exp(TA + F2){I) - (exp(rFi) exp(rF2))(/)(c), (44) 

if we set in the commutator we obtain for Fi and F2 

ErrA-BX2 = riF^ic) o Fi)(c) - Fi{c) o F2){c)) + 0{r^), (45) 

The Strang Splitting is given as: 

Errstrang = exp(rFi + ^2 ) (/) (46) 
-(exp(r/2Fi) exp(TF2) exp(T/2A))(/)(c), 

if we set in the commutator we obtain for Fi and F2 

ErrstrangA,2 = ^^'([^2, [^^2,Fi]](c) - 2[Fi, [Fi , i^^s]] (c)) + 0{t^) . (47) 
with the commutator: [Fi,F2]{c) = (F^(c) o Fi)(c) - iF{{c) o F2){c). 

5 Improving the Initialization of Operator Splitting 
Methods 

A delicate problem in splitting methods is to achieve sufficient accuracy in the 
first splitting steps, see [24]. 

Next we discuss the extension to the improvement with Zassenhaus formula. 
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5.1 Higher order A-B splitting by Initialization 

The idea is based on the novel commutator that is given in Section |4l 

Theorem 1. We solve the initial value problem by using the method given in 
equations and 113 5\) . We assume bounded and nonlinear operators A and B. 

The consistency error of the A-B splitting is 0{t), then we can improve the 
error of the A-B splitting scheme to 0(tP),p > I by improving the starting 
conditions cq as 

Co = {T^j=2 exp{Cjt^))co 

where Cj is called a nonlinear Zassenhaus exponents that is given with the novel 
commutator, e.g. the linear case is given in [?]. 

The local splitting error of A-B splitting method can be read as follows 

Pn = (exp(r„(i + B)){I) - exp(T„S)(/) exp(r„i)(/))(c:p) (48) 

where Ct is a function of nonlinear Lie brackets of A and B. 

Proof. Let us consider the subinterval [0, t], where t = t, the solution of the 
subproblem is: 

c*(t)=exp(ti)(/)(co) (49) 
after improving the initialization we have 

c*{t) = exp(ii)(/)(^J^2 exp(Qi^))(/)(co) (50) 
the solution of the subproblem (1551) becomes 

c**{t) = exp{tB){I) exp(ti)(/)(7rP^2 exp(Qt-'))(/)co (51) 
= (exp(T„(S + i))(/)(co) + 0{tP+') 
with the help of the Zassenhaus product formula. 

Remark 1. For example, the second order A-B splitting after improving the ini- 
tialization is 

c**{t) = (^exp(tB)(/) exp(ti)(/) exp(-it2[B, i])(/)) (co) 

= exp(i(B + i))(/)(co) + 0(t3) (52) 
and the third order A-B splitting after improving the initialization is 

c**{t) = (^cxp{tB){I) cxp{tA}{I) cxp{-h^B,A]){I) exp{^t^[B,[B, A]]{I) 

-i[i,[i,B]](/)))(co) (53) 

= exp{t{B + A)){I){co) + 0{t''), (54) 
where the commutator is given as [A, B]{c) — {B'{c) o A){c) — {A'{c) o B){c). 
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Remark 2. The same idea can also be done with the Strang-Sphtting method, 
see the linear case in j24j . 

6 Alternative Approaches: Iterative Schemes for 
Linearization 

In the following, we discuss the fixed point iteration and Newton's method as 
alternative approaches to linearize the nonlinear problems. 
We solve the nonlinear problem: 

F{x) = 0, (55) 

where F : R" ^ H". 



6.1 Fixed-point iteration 

The nonlinear equations can be formulated as fixed-point problems: 

X = K{x), (56) 

where K is the fixed-point map and is nonlinear, e.g. K{x) = x — F{x). 
A solution of (1571) is called fix-point of the map K. 
The fix-point iteration is given as: 

Xi+i = K{xj), (57) 

and is called nonlinear Richardson iteration, Picard iteration, or the method of 
successive substitution. 

Definition 5. Let fi < ]R" and let G : Q ^ H™. G is Lipschitz continuous on 
Q with Lipschitz constant 7 if 

|lG(a:)-G(y)|| <7lk-y||, (58) 

for all x,y ^ f2. 

For the convergence we have to assume that K he a contraction map on f] 
with Lipschitz constant 7 < 1. 

Example 1. We apply the fix-point iterative scheme to decouple the non-separable 
Hamiltonian problem. 



q» = -^(Pi,qj-i), (59) 
p« = --^(p*-i,qi), (60) 

while we have the initial condition for the fix-point iteration: 
(po,qo) = (p(i"),q(i")) 

we assume that we have convergent results after i = 1 . . . , m iterative steps 
or with the stopping criterion: 

max(||pi+i - Pill, ||qj+i - q^H) < err, 

while II • II is the Euclidean norm (or a simple vector-norm, e.g. ^2)- 
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6.2 Newton's method 

We solve the nonlinear operator equation (|55|). 

While F : D C X ^ Y with the Banach spaces X, Y is given with the norms 
II • \ \x and II • ||y. Let F be at least once continuous differentiable, further we 
assume xq is a starting solution of the unknown solution x* . 

Then the successive linearization lead to the general Newton's method: 

F'{x^)Ax, = -F{x,), (61) 

where Axi = x^+i — Xi and i = 0, 1, 2, . . . . 

The method derive the solution of a nonlinear problem by solving a sequence 
of linear problems of the same kind. 

Remark 3. The linearization methods can be applied to the MPE methods. Non- 
separable Hamiltonian problems can decoupled to separable Hamiltonians, see 
the ideas in jT4] . 



7 Numerical Examples 

In the following section, we deal with experiments to verify the benefit of our 
methods. At the beginning, we propose introductory examples to compare the 
methods. In the next examples, some ideas to applications in Burgers and Hamil- 
tonian problems, are done. 

7.1 First Example: Burgers equation 

We deal with a 2D example where we can derive an analytical solution. 

dtu = -udxU - udyU + fiidxxU + dyyu) + f{x,y,t), (62) 
{x,y,t) enx[0,T] 

u{x,y,0) = Ua,na{x,y,0), {x,y) e f2 (63) 

with u{x, y, t) = Uana(a;, y, t) ondQ x [0, T], (64) 

where Q — [0, 1] x [0, 1], T = 1.25, and /i is the viscosity. 
The analytical solution is given as 

"ana(a;, 2/, t) = (I + exp( ^"^^^ ^ ))~\ (65) 

where /(x, y, t) = 0. 

For the non-asymptotic case we compute the right-hand side as: 

) + fix, y, t) , m{x, y, t) e [0, 1] X [0, 1] x [0, 1.25] 
u{x,y,0) (x, y, 0), (x, y) G [0, 1] x [0, 1] (67) 

with y, t) — Uasym (x, y, t) ondn x [0, T] , (68) 
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We discretize with Ax, Ay = 1/40, At = 0.001. 
The operators are given as: 

A[u)u ~ —udxU — udyU, hence A{u) — —udx ~ udy (the nonlinear operator), 
Bu — ^{dxxU + dyyu) + f{x, y, t) (the linear operator). 

We apply the nonlinear Algorithm (j40|) to the first equation and obtain 

A{ui-.i)ui = -Ui^idxUi - u^^idyUi and 

BUi-l = ll{dxx + dyy)Ui^l + /, 

and we obtain linear operators, because Ui-i is known from the previous time 
step. 

In the second equation we obtain by using Algorithm (pij) : 

A{ui-i)ui = -Ui^idxUi - Ui-idyUi and 
Bui+i = ii{dxx + dyy)ui+i + /, 

and we have also linear operators. 

The maximal error at end time t — T is given as 

errniax — l^num '^ana| — Uiax | Wnum (^i , ^) ^ana(^i,^)|, 

i—1 

the numerical convergence rate is given as 

p = log(err^/2/err,i)/log(0.5). 

We have the following results, see Tables [1] for different steps in time and space 
and different viscosities, see also the results in [23] ■ 

7.2 Separable Hamiltonian 

We deal with the evolution of any dynamical variable M(q, p) (including q and 
p themselves) is given by the Poisson bracket. 

An example for a separable Hamiltonian is give as: 

Hip,q)^^ + V{q), (70) 
A and B are Lie operators, or vector fields 
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Ax = Ay 


At 


errz,! 


crrniax 
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Pmax 


1/10 


1/10 


0.0549 


0.1867 






1/20 


1/10 


0.0468 


0.1599 


0.2303 


0.2234 


1/40 


1/10 


0.0418 


0.1431 


0.1630 


0.1608 


1/10 


1/20 


0.0447 


0.1626 






1/20 


1/20 


0.0331 


0.1215 


0.4353 


0.4210 


1/40 


1/20 


0.0262 


0.0943 


0.3352 


0.3645 


1/10 


1/40 


0.0405 


0.1551 






1/20 


1/40 


0.0265 


0.1040 


0.6108 


0.5768 


1/40 


1/40 


0.0181 


0.0695 


0.5517 


0.5804 



Table 1. Numerical results for 



le Burgers equation with viscosity n = 0.05, 



initial condition Uo(0 = Cn, and two iterations per time step. 



where we have abbreviated v = p/m and a(q) = — W{c()/m. The exponential 
operators e'*'* and e^^ are then just shift operators. 

S{h) = Qh/2B^hA^h/2B 

That is also given as a Verlet-algorithm in the following scheme. 
We start with (qo,vo)* = (q(r), v(r))*: 

(qi.vi)* = e''/2^(qo,Vo)* = (/+ i/i53a(q)^)(qo,vo)* (72) 

i 

= (qo,vo + ^Mqo))*, (73) 



(q2,V2)* =e'^^(qi,vi)* = (/ + V A)(q,, vi)* (74) 

. cqi 

= (qi +/ivi,vi)*, (75) 



(q3,V3)' = e''/2S(q^^^^^t ^ (/+ l/^^„(q)^)(q2,v2)* (76) 

i 

= (q2,V2 + ^Mqi))*- (77) 

And the substitution is given the algorithm for one time-step n — > n + 1 : 

(q3, vs)* = (qo + hwo + ^o(qo), vo + ^a(qo) + ^a(qo + /ivq + ^a(qo)))*(78) 
while (q(t"+i),v(t"+i))* = (q3,V3)*. 
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Remark 4- Here, we linearize with respect to the time-steps and assume to com- 
pute a large time sequence. 

The numerical error is given with 0{h^) based on the second order approaches 
of the Strang-splitting. 

8 Conclusions and Discussions 

We have presented novel MPE approaches to nonlinear differential equations. 
Based on the ideas of iterative splitting schemes, we could linearize the numeri- 
cal scheme and apply the linear MPE approach.. Numerical examples confirm the 
applications to nonlinear equations. In the future we will focus us on the devel- 
opment of improved MPE methods with respect to non-separable Hamiltonian 
problems. 
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